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The magnetic properties of densely packed magnetic nanoparticles (MNP) assemblies are investi- 
gated from Monte Carlo simulations. The case of iron oxide nanoparticles is considered as a typical 
C/3 . example of MNP. The main focus is put on particle size and size polydispersity influences on the 

magnetization curve. The particles are modeled as uniformly magnetized spheres isolated one from 
each other by a non magnetic layer representing the organic coating. A comparison with recent 
experimental results on 7— Fe2C>3 powder samples differing by their size is given. 
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I. INTRODUCTION 



The physics and chemistry of nanoscale magnetic particles (MNP) still gives rise to an important 
research activity due both to their wide range of potential applications and their own fundamental 



interest 



-4|. Among the large variety of MNP, iron oxide based ones 7— Fe20<3 and Fe3C>4 take a 
particular place in the field of biological and medical applications because of their bio-compatibility 
and suitable superparamagnetic properties. To translate intrinsic properties of nanoparticles to 
various applications, there is a need to control nanoparticle dispersions. Consequently nanoparticles 
are usually coated by an organic surfactant in order to prevent aggregation. The influence of 
this non magnetic layer and then the nanoparticles contact distance play a major role on collective 
magnetic properties [7(. A complete understanding of the macroscopic magnetic properties of 
MNP assemblies in terms of their individual intrinsic characteristics on the one hand and of the 
size distribution and volume concentration on the other hand is of crucial importance. Indeed this 
is a mean to get informations on the relevant parameters of the distribution and MNP properties 
from the magnetic measurements. Two key features which strongly influence the macroscopic 
magnetic properties of these systems are the magnetic structure at the particle scale, where core 



8Mllj and the size distribution generally 



shell structure and spin canting effect can be invoked 
described through a lognormal law for the diameters distribution. 

At temperatures higher than the blocking temperature T& where the MNP are in the super- 
paramagnetic regime [ij, 4 and in case of weak interparticle interactions, namely for both small 
particles concentrations and in the absence of cluster formation the physical properties of MNP as- 
semblies are well understood. The magnetization curve, M(H) of the whole assembly follows then 
a Langevin like function weighted by the diameter distribution function and eventually modified in 



order to take into account a core-shell structure 



12 



14] . Moreover the one-body magnetocrystaline 



anisotropy energy of the MNP can also be taken into account and this modifies the M(H) curve 
from the original Langevin function [15|, LL6J . The core shell structure of the MNP may consist 
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ll| or of the in- 



simply of the inclusion of a magnetic dead layer at the surface of the MNP 
traduction of an additional paramagnetic component in the MNP The symmci rv breaking a! 
the surface can lead to surface effects on the anisotropy energy of each MNP with noticeable effects 
on the M(H) curve 



n 



In case of diluted assemblies of spherical MNP when the particles 



are non or weakly interacting, the non interacting particles type of approach of the magnetization 
curve leads to a reasonable determination of the characteristics of the individual particles and of 
the size distribution namely the median diameter d m and the ln(d) standard deviation a. However, 
when the NP concentration increases, the interparticles interactions must be taken into account. 
These ones which for spherical and well coated MNP include mainly the interactions between the 
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MNP magnetic dipoles (DDI), have been widely studied and a large amount of works and meth- 
ods are thus available going from mean field approximation, thermodynamic perturbation theory 
(TPT) QQ for weakly interacting systems to numerical simulations for moderate to strongly 
interacting systems [ill Isil [22 1 ■ The mean field and TPT provide an illustrative physical picture 



of the relation between the local structure and either the magnetization in terms of the applied 
field or the susceptibility. For instance the demagnetizing field effect depending on the external 

n 

shape of the system, is well reproduced by the TPT [20j . As a link between TPT and numerical 
simulations, the description based on the interaction fields distributions [23j which explains the 
DDI induced reduction of the magnetization of an isotropic system as a generalization of a similar 
result obtained using the TPT and suggests that the DDI induced reduction of the magnetization 
is not related to an antiferromagnetic behavior. However, for strongly interacting systems, as in 
lyophilized powder samples or high concentration MNP assemblies embedded in a non magnetic 
matrix the numerical simulations seem more adapted. Although numerical simulations of magnetic 
properties of MNP assemblies are now many, a systematic study of the mean size and polydispersity 
effects especially for randomly organized particles with high concentration is still missing. 

The aim of this work is to investigate this problem and to interpret recent experimental mea- 



surements 



24j on powder samples of maghemite MNP assemblies differing by their median size. 



We present a Monte Carlo (MC) simulation of the mean particle size and polydispersity effect on 
the DDI in random and densely packed spherical clusters of coated spherical maghemite MNP. 
Our main purpose is to model the case of lyophilized powders or high concentration of particles 
embedded in a non magnetic matrix. A particular attention is paid to the linear susceptibility x, 
and its dependence on the median size of the size distribution. It is found that x as a function of 
d m may present a plateau, leading to a quasi independence of the magnetization with respect to 
d m in the vicinity of the low external fields. The magneto crystalline anisotropy is then shown to 
play a role for larger values of the field when the particles remain in the superparamagnetic regime 
in agreement with the findings of Refi_ I5I \ 



201 ] and in preceding MC simulations 



magnetization curves of Ref. 



22 



25] for non interacting particles, in the TPT regime 
26j . As an application, we focus on the experimental 



24|. 



II. MODEL FOR DENSELY PACKED ASSEMBLIES 



The model we use is designed to simulate the properties of either lyophilized powders samples or 
high concentration nanoparticles assemblies embedded in non magnetic matrix. As is usually done 
to model single domain MNP, the nanoparticles are modeled as non overlapping spheres bearing 
at their center a permanent point dipole representing the uniform magnetization of the particle 
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(super spin). The moment of each particle is equal to its volume times the bulk magnetization, 
M s , which means that neither spin canting effect nor magnetic dead layer at the particle surface is 
considered. We also include the magneto crystalline anisotropy with the same anisotropy constant 
K\ on all particles. The particles are supposed to be coated by a non magnetic layer of thickness 
A/2, representing the usual coating by organic surfactant molecules. The layer thickness is taken 
as A/2 for convenience (see below). The particle diameters, {di} are distributed according to a 
log- normal law defined by the median diameter d m and the standard deviation a of ln(d), 

d m and a are related to the mean diameter and the diameter standard deviation Od through 



d\ — d m e a I 2 and a c i — d m yj (e a2 — l)e a2 . In the following, we use d m as the unit of length, 
and thus in reduced unit, the distribution function is totally determined by the single parameter 
a which characterizes the system polydispersity. 

We consider mainly spherical clusters, where owing to the global shape isotropy the demagne- 
tizing effects vanish, with free boundary conditions. This choice of large spherical clusters can be 
justified on the experimental point of view since upon drying the NP are likely to aggregate in 
spherical shaped large clusters which has been confirmed from simulations 5]. Our first purpose 
is to focus on the contribution of the dipolar interactions (DDI) to the magnetization curve, espe- 
cially in the moderate to strong coupling regime when particles surrounded by their coating layer 
are at contact. The geometrical configuration of two particles of different sizes at contact with 
their coating layer is displayed in figure ([!]). Moreover, we consider temperatures such that the 
particles of size d m are superparamagnetic; as we shall see later for polydisperse systems due to the 
presence of large particles in the distribution, this condition may not be strictly fulfilled. When 
taken into account the magnetocrystalline anisotropy is considered in its simplest form, namely 
in the uniaxial symmetry and at lowest order 2j, |4j. The total energy thus includes the DDI, the 
one-body anisotropy term and the Zeeman term corresponding to the interaction with the external 
applied field H a = H a h. Let {ri}, {v(i)} , {rhi} and {Hi} denote the particles locations, volumes, 
moments and easy axes respectively. The total energy of the cluster reads 

E= ^2^ mim ^ ^ 3 - K 1 ^v(t)(n i m i y ~ ^oHa^m^h (2) 

i<j l J i i 

where hated letters denote unit vectors, mj are the moment magnitudes, ry = \fi — fj\. It is 
worth mentioning that the consideration of the anisotropy term with a fixed easy axes distribution 
means that the magnetization relax according to a Neel process 
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28J, namely the particles are 
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considered fixed while their moment relaxes relative to their easy axis. In this work only the case 
of a random distribution of easy axes is considered. In the following we use reduced quantities; first 
the energy is written in ksTo units, To being a suitable temperature (To = 300-ftf in the present 
work) and we introduce a reference diameter, d re f. The reference diameter, d re f is a length unit 
independent of the size distribution, useful for the energy couplings, and can be chosen from a 
convenient criterion independently of the actual structure of the MNP assembly. The reduced 
total energy is given by 

fo) ( d m \ ™ ^2 (o) ( d m \ 3 \ - ,*3,j*3 ^m? - Sjrhifij^mjfij) 



r* 3 

i<3 13 



£ (0) = ^0 (7r/6) 2 M 2 d 3 e/ £ (0) = 0oKlv(dref 



dm \ Hn 



47T 

h = PoVoM s (ir/6)d 3 m H a = ( -pL ) _p- (3) 

\a re f J Ilref 

where /3q = (fe^To)" 1 and the stared lengths are in d m unit. The dimensionless dipolar coupling 
constant and anisotropy constant are then = {d m / d re f) s e^ and ck = (d m / d re f) 3 respec- 
tively; the reference diameter, d re f can be chosen such that ed{d m = dref) = e d = ^ > ^ ne re duced 
external field h coincides with the usual Langevin variable at temperature To for a monodisperse 
distribution with d = d m In equation ([3]), we also introduce the reference external field, H re f for 
convenience. 

Concerning the structure in position, the nanoparticles surrounded by their coating layer A form 
an assembly of hard spheres of effective diameters {di + A} (see figure ([T])) which are arranged in 
large densely packed clusters with either a random or a well ordered structure (simple cubic or face 
centered cubic lattice). We build these clusters in two steps. First a large stacking of the coated 
spheres is made in a parallepipedic box with the desired structure, random or well ordered. In 
the random case, this first step is made from a sequential random rain plus compression algorithm 
in such a way to maximize the packing fraction. Doing this we can get a packing fraction tp for 



the effective spheres corresponding to the so-called loose random packing 29] ((p ~ 0.60 in the 
monodisperse case). Once this first step is performed, we cut within the global stacking the cluster 
we want to study by imposing both the external shape, either spherical or prismatic, and the 
number of particles N p , with typically N p ~ 1000. The central part of some of the clusters used 
in the present work corresponding to different values of the polydispersity, a is shown in figure ([2]) 
to illustrate the structures obtained. It is important to note that because of the coating layer of 
thickness A/2 the closest distance of approach between particles i, j is shifted from (di + dj)/2 to 
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(di + dj)/2 + A and as a result the sum involved in equation @ depends on the actual magnetic 
particles concentration of the cluster through the value chosen for A. One can rewrite the DDI 
sum by using another length scale, namely (d m + A) in order to exhibit a contribution independent 
of A. Doing this, the total DDI energy reads 



(o) / d m 



3 3771,7% - iirh.fij^rhjfij) 



(rW(d ro + A))3 



(4) 



We recall that the distribution of reduced diameters, {d*} depends only on the value of a, which 
is conserved through a scaling operation corresponding to a change of d m . The sum of equation 
([3]) is a geometric sum characteristic of the DDI expected, at least for small values of a, to be 
independent of A and thus to characterize the reduced DDI sum of the most concentrated cluster 
(A = 0) of the structure (s.c, f.c.c, random) considered. In other words, equation (j4]) allows to 
explicit the dependence of the dipolar coupling with respect to the particles volume fraction, ip v . 
For this we note that (d m /(d m + A)) 3 can be rewritten as <p v /tp m where ip m (= <p v (A — 0)) is the 
maximum value of (p v for the given configuration, namely the volumic fraction corresponding to 
the spheres of diameters {di + A}. ip m = 7r/6, 0.74, and ~ 0.60 for the simple cubic, fee and the 
loose random packed structures respectively. Then from we can introduce an effective dipolar 
coupling constant, say e^^ 1 = ef\d m /d re f) 3 (d m /(d m + A)) 3 which is rewritten as ejj^ = 
{<p v / l Pm)(dm/dref) 3 , since e|jj =1.0. Now, one can replace both d m and A by say d m 2 and A2 
respectively in such a way that the total DDI energy remains constant by imposing 



(5) 



leading to 



d m 2 — 



2(1 + A/d m ) 



,A 5 



1/2' 



(6) 



In the absence of anisotropy energy, namely when only the DDI is taken into account, the two 
systems characterized by (d m , A) and (d m 2, A2) are similar and therefore present the same mag- 
netization curve in terms of the reduced field h. Furthermore this holds also whatever the value 
of ck in the vicinity of zero external field because for random distribution of easy axes the linear 
susceptibility x does not depend on ex in the superparamagnetic regime. Doing the transformation 
©, the actual values of {fi} are scaled according to the value of (d m + A). Our hypothesis of 
a value of a for the reduced diameter distribution to be not (or only negligibly) modified holds 
rigorously in the quasi monodisperse case (a << 1). Consequently we shall use in the following 
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the scaling transformation © only in quasi monodisperse situations. 

In the present work we focus on both the reduced magnetization per unit magnetic volume in 
the direction of the external applied field, 



Ms £i"(0 



and the linear susceptibility, 



8M M s fd m \ 3 dM r _ M s fd m \ 3 8M r (h) 
X -M-H^{-d-])^h-H^{d^) Xr Wlth Xr ~^ir (8) 

where we have used equation ^ to introduce the reduced susceptibility, Xr- The susceptibility 
can also be obtained from the fluctuations : 

^ v(i) ( < > , 

xr = ^E-7rr \v cJ ' (9) 

*-r J v[d m ) I (Ei«W) ' 

As a rule, we use this second way with the direct derivative merely used as a check of the calculation. 

When the anisotropy energy is zero, the magnetization curve can be simulated either starting 
from h = and increasing the field step by step or from the starurated situation, and decreasing 
h down to h = 0. When the anisotropy energy is included and since we may get an opening of 
the hysteresis loop, we start from the saturated case at sufficiently high applied field, and decrease 
the field beyond —hi rr where the irreversible field hi rr is defined as the value of h below which the 
hysteresis cycle opens. In cases where the hysteresis cycle opens, we also define an anhysteretic 
magnetization curve from the downward and the upward magnetization curves which because of 
the symmetry of our system reads 

M^ an \h) = i(M r W(/i) + MW(li)) = i(M r W(/i) - M«(4)) (10) 

The magnetization curves M(h) in terms of the external field are determined from Monte Carlo 
simulations, by fixing the locations of the particles in the cluster. We consider free boundary con- 
ditions, and the clusters includes ca 1000 particles. The dipolar coupling parameter is determined 
from equation @- In section fill CI we consider a given set of experimental results in order to illus- 
trate the model; nevertheless we do not restrict this latter only to this well specified set of samples 
but instead use the characteristics of maghemite as typical example for MNP assemblies. For the 
bulk magnetization M s we use a commonly accepted value for maghemite. Using M s varying from 
80 to 84 emu/g, or ~ 75.0 emu/g, if we take into account the temperature dependence, and p 
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= 4.870g/cm 3 leads to fi M s from 0.459 T to 0.514 T; at T = 300K we get 4° } = 1.0 for d ref 
varying from 9.665 nm to 10.422 nm and we use in the following except otherwise mentioned d re f 
= 10 nm which corresponds to fioMs = 0.488T and H re f = 16.20 kA/m. The anisotropy constant 
K\ cannot be taken equal to the bulk effective magnetocrystalline anisotropy constant Kb as it is 
found to be much larger when the particle size decreases. A rather wide spectrum of values can be 
found in the literature for K±, corres pon ding to (K\jKb) lying in between ~ 4 to 15 for particle 



diameters of ca 12 nm or smaller 



30j-[34|. In the following we use either (Ki/Kb) ~ 4 or 2, since 



we consider particles with mean diameters larger than 10 nm. With Kb = 0.47 10 4 J/m 3 [10l.l3ll 
this leads to — 2.38 and 1.19 respectively. In any case, both M s and K\ are to be understood 
merely as realistic instead of truly accurate experimental values given the simplicity of the model. 
Our Monte Carlo simulations are based on the usual Metropolis scheme 35j, [36( ; the averages 



are taken over 10 to 40 independent runs each of which consists in 10 4 to 2.10 4 MC steps (MCS) of 
equilibration followed by 2.10 4 to 3.10 4 MCS for the averages calculations. Each MCS consists in 
one trial move per moment in average. The trial move on the unit moment rhi consists in moving 
rhi to (rhi + wu) / \rhi + wu\ where u is a random vector picked within the unit sphere with uniform 
probability density. This remains to move rhi in a cone of maximum deviation 89 whose value is 
controlled by the amplitude parameter, w. For 89 << 1, we have 89 ~ w. The value of w 
can be either fixed for a time scale mapping of the MCS or determined in a self consistent way 
in order to optimize the sampling by imposing a value for the acceptance ratio, R. The former 
version of this scheme corresponds to the time quantified Monte Carlo algorithm (TQMC) 37, 3|| 
in its first formulation ignoring the precessional step 37j . In the absence of anisotropy energy, the 
time scale mapping is irrelevant for the present purpose since we expect neither a ferromagnetic 
behavior nor a metastable blocked regime. Thus in this case, w is self consistently determined in 
such a way that R — 0.5. Conversely, when e^ 1 ^ 0, especially for polydisperse distributions 
we expect the largest particles to be in blocked state leading to a remanent state all the more that 
the DDI increase the blocking temperature. Hence, especially in the vicinity of h = 0, we deal 
with a metastable state whose life time must be comparable to the long scale measuring time r m . 
Strictly speaking one has to perform MC simulations corresponding to r m and to use the version 
of the scheme outlined above allowing a mapping of the MC step on the true relaxing time. Since 
we are interested only in the long time behavior (corresponding to the SQUID measurements time 
scale), we do not focus on a precise mapping of the MCS scaling time. Instead, we determine w 
from the behavior of the instantaneous polarization M(t), versus t in MCS along a MC run at 
h = starting from {rhi} = z. In other words, we chose w in order to avoid nonphysical jumps 
over the anisotropy energy barrier. By varying w we get as expected a w dependent evolution of 
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M(t) before reaching a fluctuating behavior around a well defined plateau; the long time mean 
value < M(t) > determined beyond some threshold t value and for t up to 2.f0 5 MCS is found 
independent of w at least for w varying in the range w — 0.03 to 0.25 for typical values of the 
parameters we consider (e^ ~ 2 to 8, ej? ~ 2.3) and the polydispersity deviation a — 0.28. 
Therefore, in the following, we Gxw = 0.25 when ^ 0. 



III. RESULTS AND DISCUSSION 



A. Weak coupling case 



Before focusing on the powder like situation characterized by a moderate to strong dipolar 
coupling, we consider the weak coupling limit of the DDI, < 1 with €k = where one can 
compare the results to the analytical one obtained from the thermodynamic perturbation theory 
and make the link with the mean field approximation. The important point is that one can 
deduce at least qualitatively when ed deviates from the limit td << 1, the general behavior of the 
magnetization with respect to the DDI. In this framework, we can expand both the magnetization 
M r (h) and the susceptibility \ r in terms of €d Qj- 

M r (h) = M^(h) + M^(h)e d + ^M^e 2 d 

Xr = X^ + X^ed+^el (11) 

Mr°' and Xr°^ correspond to the non interacting case, namely 



(0) _ Jv(d)L((3*(d/d m fh)f(d)d(d) 
r W ~ fv(d)f(d)d(d) 



(12) 



where L is the Langevin function; Xr°^ is directly related to Mr°\h) and at h = leads to the 
linear susceptibility of the non interacting system 

X < ?H0)=P*(dt/d* 3 )/3 = /Texp(27a 2 /2)/3 (13) 



where d* is the s-th reduced moment of the distribution f(d). Equation (fT3j) explicits the effect of 
the polydispersity through the factor (dg/cfij), written here in terms of a for the lognormal distri- 
bution. The expansions ([lip which have been explicited in the framework of the TPT in 
depend on geometrical sums which can be directly calculated from the structure considered. More- 
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over, the linear part with respect to h of Mr \h) can be deduced in the mean field approximation 
of the magnetization which introduces the DDI contribution to M r from the demagnetizing field 
and follows from M r {H a ) — Mr°\H e ft) where Mr°' corresponds to the non interacting system 
and H e f f is the effective field 

H eff =H a - (D a - l/p(°)(ff a ) ; (14) 

D a is the demagnetizing factor of the sample in the direction of the external field, h — e a , and M v 
is the total magnetization per unit volume which is related to either the number of MNP per unit 
volume, p, or the MNP volumic fraction, <p v , through 

= M s M^(H a )p J v(d)f(d) = M s M^(H a )^ v (15) 

Using equation (Tl2|) for M,^ ' and keeping only the first order term with respect to h, we get 

H eff =H a -\(D a -l/3)M s l3*^h (16) 

which is then inserted in the mean field expression for M r (H a ); then form an expansion of M< 0) 
at first order with respect to h and from equation © for ed, we get 

2 



MV{h) = -f)*\D a - 1/3)^Y|) h (17) 



Equation (|17[) can be equivalently rewritten, in terms of AM r = M r {td) — M r (ed = 0), as 



d(AM r (h)) _ „ a ._ i/a .8 (- fd 6 



2 



*ndC 2 = - ' {h = 0,e d = 0) = ^- = -P'(D a -l/3)-<p v [-± 



(18) 



A result in agreement with refs. 



Q and 3 



in the monodisperse case. Here, the important point 



is that we explicitly write down the effect of the polydispersity through the factor (d§/d%) 2 which 
strongly deviates from unity once a takes non negligible values. It is worth mentioning that the 
preceding equations hold when either e<j << 1 or (p v << 1. We have performed MC simulations 
of the magnetization at small values of the coupling constant for prismatic clusters corresponding 
to either well ordered (simple cubic, and c.f.c) or random structures with a monodisperse particles 
distribution, and a random structure with a polydispersity characterized by a — 0.28. The results 
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for the second derivative of AM r with respect to td and h, C2, is displayed in table U As can be 
seen, especially for h = z, the mean field approximation or equivalently the linear contribution 
of the TPT compares well with the MC simulations and in particular for the polydispersity effect. 

Since spherical or cubic systems are characterized by D a = 1/3, this first term vanishes in these 
situations and one is left for the DDI contribution with AM r oc e 2 d and similarly for Xr- Moreover, 
still for isotropic systems, we know that the DDI contribution to both M r and x is negative. 
Therefore the magnetization is all the more reduced due to the DDI that the coupling constant €d 
increases. From the analytical results of the TPT we can calculate the proportionality coefficient 
relating (xr(td) — Xr(0)) to A. We have thus compared the MC simulation to the theoretical 
small ed expansion in the simple cubic structure case and a monodisperse distribution. From this 
comparison, see figure EH we can check that the TPT gives an accurate result only for < 0.2 
as expected. Furthermore, from a description based on the dipolar fields distributions which can 
be seen as a generalization of the mean field type of approach, ref. [23J have shown also that 



the dipolar interactions in isotropic systems decrease the magnetization. This decreases is related 
to the non linearity with respect to the applied field of the non interacting contribution Xr°^ to 
the susceptibility. Notice that this second type of approach, which remains qualitative in the 
absence of a theory to deduce the dipolar field distribution, is not restricted to the weak coupling 
case. Hence, as a general rule, we expect that in an isotropic sample the DDI tend to reduce 
the magnetization. However, this reasoning does not hold at high fields where the Zeeman term 
dominates on the DDI and where we expect an approach to saturation, close to what is obtained 
in the non interacting case deduced from the high field expansion of M^ ' of equation (|12p , namely 
(M r (h) ~ 1- \/{p*dlh)). 



B. Spherical clusters in the strong coupling case 

We now consider, exclusively for spherical clusters, the moderate to strong dipolar couplin. 



24] 



case corresponding to the experimental situation of typical coated maghemite NP powders 
with e d = 1 for d re f = IQnm and a coating layer of A = 2.0nm. The median diameter varies 
from d re f to 2 x d re f and the standard deviation of the distribution ln(d) is taken from a = 0.05 
to represent the quasi monodisperse case to a — 0.50 to represent a large polydispersity. The 
importance of a on the MNP distribution in the clusters is clearly seen on figure Notice that 
a standard value obtained experimentally is ca. 0.20 ~ 0.30 which is represented here by a = 
0.28. In the first step we neglect the anisotropy contribution (ck = 0) and focus only on the DDI. 
First of all we analyse the linear susceptibility, x which provides the behavior at low field of the 
magnetization. Since in our model, with a constant coating layer thickness, A, the dipolar coupling 
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constant scales as (d m /d re f) 3 we expect in the vicinity of h ~ a reduction of the magnetization 
higher for large median diameters, where the initial non interacting magnetization M^ ' is higher. 
In the quasi monodisperse case, a = 0.05 we make use of the scaling transformation introduced 
in equation ([5]) to explicit the effect of the coating layer thickness A on x by using only one set 
of simulations for A/d re f = 0.20. We checked for A = 0.8 and 2 values of d m the reliability 
of this scaling transformation (see figure (JIJ)). Therefore, in the quasi monodisperse case we have 
a rather complete picture of both the effect of the variation of the median diameter, d m and of 
the distance of closest approach between NP, controlled by the coating layer thickness, A. The 
result for \ determined with a = 0.05 is displayed and compared to the non interacting case x*- -* 
in figure |4]). As expected, when increasing the particle size and consequently the DDI coupling 
constant e<i, an increasing reduction of x is obtained. This reduction is of course all the more 
important that A is small. The important result is that we can get a plateau, which means that 
X may becomes particle size independent beyond a threshold value which is, as expected, strongly 
A dependent. As already mentioned, for random distribution of easy axes, x does not depend on 
ex in the superparamagnetic regime, and accordingly this result holds also in the case where the 
anisotropy is included. 

The dependence of x on both d m and A can be used to deduce the behavior of x with the NP 
volumic fraction (or concentration) at fixed value of d m through the relation ip v = y> m ( d m /(d m + 



22, 



39(, we get a 



A)) 3 with (p m — (p v (A = 0). Doing this, in agreement with other MC results 
monotonous decrease of x with the increase in ip, as shown on figure ([5]). Furthermore this shows 
that a fit of the NP size on the magnetization curve by using a Langevin function does not hold 
beyond a critical value ip c of the volumic fraction. We can estimate this latter by imposing that 
x/x( e d — 0) is larger than some threshold value say A, leading the determination of tp c through 
x{Vc)/x{ e d = 0) = A. The result obtained by using A = 0.80 is displayed on figure 

The magnetization curves in terms of the reduced external field H a / H re f for three values of the 
median diameter, still for a = 0.05 is shown in figure © and compared to the non interacting 
diameter distribution weighted Langevin curves. We clearly see the important reduction of the 
magnetization compared to the non interacting case, and the very weak dependence of the low 
field behavior with respect to the median diameter which is expected as the considered sizes are 
either close to the onset of the x(dm) curve plateau corresponding to (A/d re f) — 0.20 [d m /d re f — 
1) or pertain to this later (d m /d re f — 1.33 and 2.00). On the other hand at low external fields the 
nearly size independence of the magnetization is correlated with a quasi linear behavior of M r {h) 
with respect to h, which seems coherent with the interaction fields distribution description |23l ] . 

Then we introduce the polydispersity at fixed values of d m . First we consider the case d m /d re f 
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= 1.33, as an example of median diameter located in the plateau region of the x(d m ) curve. In this 
case we expect a very weak dependence of the magnetization with respect to the polydispersity 
in the low field region and this is confirmed by the MC simulations. Indeed, we get only small 
changes of M r (h) with a as can be seen in figure ([7|). The magnetization curves corresponding to a 
up to 0.40 are very close to each other for the values of the field for which M r < 0.70; beyond this 
value, the deviations between the different magnetization curves reflect mainly the approach to 
saturation where M r (h) ~ (1 — l/(/3*d^h)) depends on a through d%. The deviation from the quasi 
monodisperse situation over the whole field range occurs for a > 0.5. Conversely, when the median 
diameter is taken outside of the xiAm) plateau, as is the case for d m /d re f = 1.0 the polydispersity 
has a noticeable influence on the magnetization as shown in figure ([8]) for a ranging from 0.05 to 
0.40. 

In the superparamagnetic regime the MNP anisotropy energy modifies the magnetization curve 
for intermediate values of the field and leads to a reduction of M r since the moments tend to 
be pinned in the easy axes directions. Taking into account €k thus reduces further M r for h 
between the low field region controlled by the DDI and the approach to saturation controlled 
by the Zeeman energy. In the quasi monodisperse case, the blocking temperature corresponds 
to that of the median diameter, namely for non interacting particles, ksTb ~ Kiv(d m )/25 or 
equivalently for the reduced blocking temperature l//3£ — (d m /d re f) 3 /25 leading to l//3£ — 
0.225 for d m /d re f — 1.33. Here we restrict to the room temperature, /3* — 1, and we expect the 
system to be in the superparamagnetic regime even for short times. Indeed for er = 0.05 our MC 
simulations confirm the superparamagnetic regime. The result is displayed and compared to the 
€k = case in figure © for d m /d re f = 1.0 and 1.33. As expected, the anisotropy energy does 
not affect the M(H) curve in the vicinity of H = due to the random distribution of easy axes. 
Moreover, when d m /d re f — 1, the M(H) curve for intermediate values of the field is only weakly 
modified by the anisotropy energy while for d m /d re f = 1.33 a noticeable deviation is obtained. 

The influence of the polydispersity on the magnetization curve when the anisotropy energy is 
included is shown for d m /d re f = 1.33 on figure (fTU|) for a ranging from 0.05 to 0.35. Because of 
the largest particles in the distribution, the system is no more in the superparamagnetic regime 
for the MC runs considered up to 10 5 MC steps. On the qualitative point of view this is expected 
since 1/ (31 behaves as d* 3 in the absence of DDI and moreover increases with the DDI. As a result 
an opening of the hysteresis cycle is obtained with remanence magnetization and coercive field 
increasing with a as shown in figure (|10p in the particular case d m /d re f = 1.33. The magnitude of 
the hysteresis cycle opening is expected to increase with d m /d re f and is indeed found very weak 
for d m /d re f = 1.0. The determination of the remanence in terms of d m and the measuring time 
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is beyond the scope of this work; we nevertheless note (see section Ull C[) that the hysteresis cycle 
opening for large values of d m is in qualitative agreement with experiment. 



C. Comparison with experiment 



We now consider experimental results obtained recently on 7 — Fe20% NP powders samples 



differing by their size 



2J. The 



24j . The experimental protocol for the synthesis is described in 
particles are coated with (5-hydroxy-5,5-bis(phosphono)pentanoic acid) which provides a coating 
layer of thickness ca 2 nm between particles. As a result of the synthesis method, the standard 
deviation of the diameter distributions as determined by TEM takes nearly the same value in the 4 
samples considered, namely, a ~ 0.26. The saturation magnetization M s is found to be in between 
61 and 70 emu/g for the distributions characterized by d m = 10, 12, 18 and 21 nm. Although 
these values are smaller than the bulk value at room temperature (~ 75 emu/g) the difference is 
small enough for the spin canting to be neglected in first approximation. It is worth mentioning 
that the magnetic properties of these NP assemblies have been measured also in diluted solution 
and, although the possible formation of clusters and/or chains in the presence of the external field 
cannot be ruled out, this allows for an estimation of the interaction effect. When going from the 
dispersed samples to the powder ones, we observe both a strong reduction of the magnetization 
and its weak size dependence in the low field region 24|. According to our simulations, both effects 



result from the DDI. The DDI induced reduction of M r (H) in the absence of demagnetizing effects 
is a general simulation result 
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39 



40j and a similar trend has been obtained experimentally 



41| 



and can be deduced from the FC/ZFC measurements in the superparamagnetic regime of either 



bare or Si coated 7— Fe2C>3 NP [33j. Beside its rather weak size dependence the other feature of 
the experimental reduced magnetization curves, M r (H) in the low field region (see figure (|11|) is 
the opening of the hysteresis cycle for the largest sizes beyond d m — 12 nm. These two points are 
in qualitative agreement with the MC simulations on our model although the opening of the cycle 
becomes noticeable for larger median diameters (d m ~ 18 nm) than in MC simulations. 

In the present work, we do not compare the experimental magnetization curve in the whole range 
of field with the results of either a mean field approach or the TPT. In any case the values of the 
dipolar coupling corresponding to the experimental powders samples (e^^ ~ 0.6 to 6.0 when 
ed ~ 1 to 8 and the effect of A is taken into account) fall outside of the range of validity of the 



TPT. Indeed this later is limited to ca < 1/6 according to ref. [20j, the analytical approach 
based on TPT of ref. [18j is shown to be very accurate for < 0.25 and valid for < 0.50 in 
the monodisperse case and in section Ull Al we found that \ as calculated from the TPT start to 
deviate from the simulated results at ed ~ 0.2. Moreover the accounting of the polydispersity is 
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expected to worsen the lack of accuracy of the TPT with the increase of 6d- 

In figure (|12[) we compare the experimental and simulated M r (H) for the applied field in the low 
to intermediate range for d m = 10 nm. The agreement is quite satisfactory up to H = 60 kA/m 
where M r ~ 0.7. Then for median diameters d m > 12 nm (d m /d re f > 1.2, we get for both 
the experimental samples and the MC simulations an opening of the hysteresis cycle. However, as 
can be deduced from figures (fTT|) and (fTU|) the irreversible field, Hi rr is found much larger in the 
MC simulations than in the experiment. Notice that we do not try to map rigorously the MC time 
scale to the actual measurement time r m , and this is plays a central role on this point. Therefore 
concerning the MC simulations we consider in the following the anhysteretic magnetization as 
defined in equation (|10j) ; using this later remains to ignore the hysteresis cycle (i.e. the remanencc 
and the coercive field) or to consider the infinite time scale limit. We compare the experimental 
M(H) to the simulated ones in figures (flUl [T4"|) for the median sizes d m 12 and 20 nm respectively. 
Given the simplicity of the model which does not include at all any structure in the particles, and 
the absence of fitting parameter the agreement is qualitatively satisfactory when = 2.38 in 
particular for the overall variation of M r (H) at low fields (H < 40 kA/m). 

Since in this range of fields, the deviation of M r from the non interacting case is governed by 
the DDI, we can conclude that the strong reduction of the M r (H) variation with respect to H 
and its relative size independence when compared to the diluted solutions counterpart is the DDI 
signature. For median sizes larger than d m = 10 nm, the main discrepancy between the simulated 
and experimental M r (H) curves is the strong non linearity in the very vicinity of H — 0. This is 
clearly due to the oversimplification of the model in which the particles are uniform single domain 
ones. 

D. Conclusion 

In this work we have used MC simulations to investigate both the median size and polydispersity 
effects on the magnetization curve of densely packed clusters of single domain magnetic NP. An 
important result is the plateau in the x(dm) curve in the quasi monodisperse case for small values 
of the coating layer A, which emphasizes the much reduced size dependence of the M(H) low field 
dependence in the concentrated systems. Despite of the simplicity of the model, some important 
features of the experimental M(H) on powder samples are reproduced, especially concerning the 
DDI signature which occurs principally at low fields and its dependence on the particle size. In 
order to get a satisfactory agreement with experiments, it appears that the internal polarization 
structure of the NP should be introduced. 
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Tables 



struct. 


cs^ 


cfcW 


rand, a = {c) 


rand, a = 0.28^ 




tt/6 


0.74 


0.525 


0.580 


& 


0.2775 


0.3367 


0.2534 


1.7934 




-0.4745 


-0.6109 


-0.4390 


-3.1506 


dp 


0.2217 


0.3133 


0.2217 


1.6173 


cf 


-0.4433 


-0.6265 


-0.4433 


-3.2354 



Table I: Second derivative C2 = d 2 (AM r ) / deddh at €d = and h = for a primatic cluster of 1024 < - a - ) , 
1103 (6) , 1054 (c) and 879 (d) particles, characterized by L x = L y and L z = L x /5. MC simulation with 
h = x (/) , h = z ' 9) ; equation (jTHJ) with h — x * h) , h — z w . w The demagnetizing factor entering 
equation (1181) is taken from [42| and the moments d* are taken from the actual diameters distribution of 
the cluster considered. 
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Figure captions 

Figure [T] Shematic view of the configuration for two particles coated by the layer of thickness A/2 at 
contact. 

Figure [2] Central part of the clusters corresponding to d m /d re f = 1.33 and a = 0.05, 0.28 and 0.50 
from top to bottom. 

Figure [3] Comparison of \ r in terms of as calculated from the TPT of Ref. [l5[ , (solid line) and the 
present MC simulation (symbols) for a spherical cluster of simple cubic structure with N p = 
1021 particles and a monodisperse distribution. (3* — 1.0. 

Figure U Linear susceptibility in terms of the median size d m /d re f for different values of the coating 
layer thickness A in the quasi monodisperse case, a — 0.05 and (3* = 1.0. The two crosses 
on the A/d re f = 0.8 curve correspond to the direct calculation without using the scaling 
transformation @. The dotted lines are guides to the eye and the solid line corresponds to 
the non interacting case. 

Figure [5] Susceptibility x in terms of the reduced volumic fraction, for a = 0.05, (3* — 1.0, and d m /d re f 
= 2.0 (solid circles); 1.50 (solid squares); 1.25 (upward triangles) and 1.0 (downward trian- 
gles). (p m = ip(A = 0). In the present work <p m ~ 0.585. Inset: reduced critical volumic 
fraction defined as x( < / 5 j e d)/x( e d = 0) = 0.80 in terms of the median particle size. 

Figure H] Magnetization in terms of the applied field for different values of the median diameter, A / d re f 
= 0.20, j3* — 1.0 and N p = 1007 in the quasi monodisperse case, a = 0.05. The corresponding 
non interacting curves (diameter distribution weighted Langevin curves) for d m /d re f = 1.0 
(long dash), 1.33 (short dash) and 2.0 (dotted line) are displayed for comparison. 

Figure [7] Magnetization in terms of the applied field for different values of the ln(d) standard deviation 
a and d m /d ref = 1.33, A/d ref = 0.20 and (3* = 1.0. N p = 1007 (a = 0.05), 923 (a = 0.28); 
985 (<7 = 0.40) and 990 (a = 0.50). The dotted lines are quides to the eye. The solid lines 
are the asymptotic limits for a = 0.05 (bottom) and a = 0.50 (top). 

Figure H] Magnetization in terms of the applied field for a — 0.05, 0.20, 0.28 and 0.40 from bottom to 
top. and d m /d re f — 1.00, A/d re j = 0.20, (3* — 1.0. The dotted lines are guides to the eyes 
and the solid lines are the corresponding asymptotic limits for a = 0.05 and a — 0.40. 

Figure |9] Magnetization curve with the anisotropy energy = 2.38, a = 0.05 and the particles sizes 
as indicated. The solid lines correspond to ex = 0. f3* = 1 and N p = 1007. 

Figure ITU1 Magnetization curve with the anisotropy energy = 2.38, d m /d re f = 1.33 and /3* = 1 in the 
polydisperse case with a = 0.35 (open circles); 0.28 (triangles); 0.20 (squares) compared to the 
quasi monodisperse case, a = 0.05 (solid line). Insert: detail of the downward magnetization 
curve in the vicinity of h — 0, showing the evolution of the remanence and coercivity with a 
for a = 0.35, 0.28, 0.20 and 0.10 (open squares) from top to bottom. 

Figure [TT] Experimental low field magnetization curves at room temperature for powder samples of 
7— FS2O3 from ref. [2~4| with polydispersity a ~ 0.26 and median sizes as indicated. 

Figure [12] Comparison between MC simulations and the experimental magnetization curve at room 
temperature for d m — 10 nm. The simulation are performed with j3* — 1, = 2.38 and 
N p = 923. 

Figure [T3] Comparison between MC simulations and the experimental magnetization curve at room 
temperature for d m = 12 nm. The simulations are performed with (3* — 1, N p — 980, ej^ = 
2.38 (open circles) or 1.19 (dashed line). The dooted line is a guide to the eye. 

Figure [T^] Comparison between MC simulations and the experimental magnetization curve at room 
temperature for d m — 21 nm. The simulations are performed with d m = 20 nm, f3* = 1, N p 

= 998, = 2.38 (open circles) or 1.19 (dashed line). 




Figure 2: Central part of the clusters corresponding to d m /d re f — 1.33 and a — 0.05, 0.28 and 0.50 from 
top to bottom. 
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Figure 3: Comparison of \r in terms of as calculated from the TPT of Ref. (solid line) and the 

present MC simulation (symbols) for a spherical cluster of simple cubic structure with N p = 1021 particles 
and a monodisperse distribution. /3* — 1.0. 
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Figure 4: Linear susceptibility in terms of the median size d m /d re f for different values of the coating layer 
thickness A in the quasi monodisperse case, a = 0.05 and /3* = 1.0. The two crosses on the A/d re f — 0.8 
curve correspond to the direct calculation without using the scaling transformation ©. The dotted lines 
are guides to the eye and the solid line corresponds to the non interacting case. 
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Figure 5: Susceptibility \ in terms of the reduced volumic fraction, for a = 0.05, j3* = 1.0, and d m /d re f 
= 2.0 (solid circles); 1.50 (solid squares); 1.25 (upward triangles) and 1.0 (downward triangles). (p m 
= ip(A — 0). In the present work ip m ~ 0.59. Inset: reduced critical volumic fraction defined as 
X(<P, £d)/x( e d = 0) = 0.80 in terms of the median particle size. 




Figure 6: Magnetization in terms of the applied field for different values of the median diameter, A/d re f = 
0.20, P* — 1.0 and N p = 1007 in the quasi monodisperse case, a = 0.05. The corresponding non interacting 
curves (diameter distribution weighted Langevin curves) for d m /d re f = 1.0 (long dash), 1.33 (short dash) 
and 2.0 (dotted line) are displayed for comparison. 
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Figure 7: Magnetization in terms of the applied field for different values of the ln(d) standard deviation 
a and d m /d reJ = 1.33, A/d ref = 0.20 and /?* = 1.0. iV p = 1007 (a = 0.05), 923 (a = 0.28); 985 (a = 
0.40) and 990 (a = 0.50). The dotted lines are quides to the eye. The solid lines are the asymptotic limits 
for a — 0.05 (bottom) and a — 0.50 (top). 




Figure 8: Magnetization in terms of the applied field for a = 0.05, 0.20, 0.28 and 0.40 from bottom to 
top. and d m /dref = 1.00, A/d re f = 0.20, fi* = 1.0. The dotted lines are guides to the eyes and the solid 
lines are the corresponding asymptotic limits for a — 0.05 and a — 0.40. 
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Figure 10: Magnetization curve with the anisotropy energy e K = 2.38, d m /d re f — 1.33 and /?* = 1 in 
the polydisperse case with a = 0.35 (open circles); 0.28 (triangles); 0.20 (squares) compared to the quasi 
monodisperse case, a = 0.05 (solid line). Insert: detail of the downward magnetization curve in the vicinity 
of h = 0, showing the evolution of the remanence and coercivity with a for a = 0.35, 0.28, 0.20 and 0.10 
(open squares) from top to bottom. 
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Figure 11: Experimental low field magnetization curves at room temperature for powder samples of 
7— Fe203 from ref. [24[ with polydispersity a ~ 0.26 and median sizes as indicated. 
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Figure 12: Comparison between MC simulations and the experimental magnetization curve at room 
temperature for d m = 10 nm. The simulation are performed with /3* = 1, = 2.38 and N p = 923. 
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Figure 13: Comparison between MC simulations and the experimental magnetization curve at room 
temperature for d m = 12 nm. The simulations are performed with fi* = 1, N p = 980, e^' = 2.38 (open 
circles) or 1.19 (dashed line). The dooted line is a guide to the eye. 
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Figure 14: Comparison between MC simulations and the experimental magnetization curve at room 
temperature for d m — 21 nm. The simulations are performed with d m — 20 nm, ft* = 1, N p = 998, e^ 1 ' 
= 2.38 (open circles) or 1.19 (dashed line). 



